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Abstract 

In this paper, we present a calculation of seismic scattering from irregular surface topography in 
the Born approximation. Based on US-wide topographic data, we investigate topographic scatter- 
ing at specific sites to demonstrate its impact on Newtonian-noise estimation and subtraction for 
future gravitational-wave detectors. We find that topographic scattering at a comparatively flat 
site in Oregon would not pose any problems, whereas scattering at a second site in Montana leads to 
significant broadening of wave amplitudes in wavenumber space that would make Newtonian-noise 
subtraction very challenging. Therefore, it is shown that topographic scattering should be included 
as criterion in the site-selection process of future low-frequency gravitational-wave detectors. 



I. INTRODUCTION 



Seismic disturbances are an important limitation to the sensitivity of ground-based 
gravitational-wave (GW) detectors. Sophisticated seismic isolation systems will be used 
in next-generation detectors like Advanced LIGO [Ij, Advanced Virgo [2j or LCGT [3j to 
push the sensitive band down to frequencies around 10 Hz. The goal of future-generation 
detectors like the European Einstein Telescope [U [5] will be to extend the detection band 
to even lower frequencies like 3 Hz and below. The detectors need to be stabilized against 
occasional seismic disturbances, as for example those generated by earthquakes, but the am- 
bient seismic field also generates a stationary background noise that cannot be fully avoided 
despite multi-stage passive and active isolation systems. 

Seismic fields also generate perturbations of the gravity field that will contribute to the 
instrumental noise as Newtonian noise [6-9j. These perturbations are extremely weak, and 
they will only be significant noise contributions in future detectors beyond Advanced LIGO 
and Advanced Virgo ^Oj . Whereas it is possible to conceive further improvements to current 
seismic isolation designs, it is not possible to shield the motion of the suspended test masses 
from fiuct nations of the gravity field. The best strategy is to search for detector sites 
that have a weak seismic spectrum within the relevant frequency band between 1 Hz and 
20 Hz [llj. In addition, because the sources of the gravity perturbations are known, it 
should be possible to estimate Newtonian noise using data from seismic arrays and then 
subtract the estimate coherently from the GW data [121 IE] • Seismic fields are difficult to 
analyze at the surface due to heterogeneities of the near-surface layers [14j and abundance 
of seismic sources. Therefore, constructing future detectors underground would not only 
provide remote, low-noise environments for the detectors but also improve and simplify the 
Newtonian-noise estimation, potentially yielding smaller subtraction residuals. 

In this paper, we address a specific problem that is related to the Newtonian-noise esti- 
mation. The claimed advantage of underground sites is partially based on the assumption 
that the medium is more homogeneous at greater depths, and local seismic sources would 
be rare or weak. This is certainly the case if underground sites are chosen carefully, but 
another effect that causes heterogeneity in the seismic field is the scattering of seismic fields 
from an irregular surface topography. We will first outline the theory of seismic scattering in 
the Born approximation in section [TT] and explicitly evaluate scattering coefficients for a few 
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examples. We then present topographical data for the US in section |TTT] and use a simple 
figure-of-merit to select interesting sites. A detailed scattering analysis is carried out for 



these sites in section [IV| The paper concludes in section |V[ where we discuss implications 
of our results for Newtonian-noise estimation and subtraction. 



II. TOPOGRAPHIC SCATTERING: THEORY 

The theory of seismic-wave topographic scattering in the Born approximation was first 
developed for the 2D case in [15j and more specifically for sinusoidal surfaces in [16]. Later, 
the 3D case was discussed in [17j. The Born approximation is convenient for the purpose of 
this paper, but it should be noted that alternative approximations have been developed that 
hold for a wider class of topographies (see [18j for details). In this section, we will outline the 
basic ideas and fundamental equations. We will consider sinusoidal surfaces in some detail, 
applying a more efficient approach compared to [16j by formulating the problem in terms 
of boundary conditions instead of effective source distributions at the surface. Scattering 
coefficients will then be plotted for the most common types of incident plane waves. 

The scenario that we consider is that of plane seismic waves propagating on or inside a 
homogeneous medium with arbitrary surface topography. If the surface was ffat, then surface 
waves would propagate without conversion into other seismic modes. Seismic waves propa- 
gating through the medium and reffecting from a ffat surface experience partial conversion 
between compressional and shear modes. In each case, we denote the total displacement 
field throughout the entire medium by ^^{x^t). The superscript "0" indicates that we con- 
sider the displacement field assuming a fiat surface. The waves emanating from a plane 
surface are called specular refiection. As we will show, an incident field is scattered from 
an irregular surface topography. In this case, we decompose the total displacement field 
^ into its unperturbed component and the scattered field We assume that further 
scattering of a scattered wave is negligible, i. e. the calculations are carried out in the Born 
approximation. 

The surface topography shall be described by a function z = s{x^y)^ which is treated as a 
perturbation of a fiat surface such that the mean value vanishes: (^(x, y)) = 0. The traction 
at a free surface, which is the stress tensor projected onto the normal direction n = riiCi of 
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the surface, vanishes (see for example [E]): 

TijUi = Tnj = (1) 

We desire to translate these boundary conditions defined on a rough surface to equivalent 
conditions on the (virtual) plane z = 0. The new boundary conditions must involve non- 
vanishing surface traction depending on s{x^y). To calculate the surface traction on the 
plane, we first rotate the traction of the free surface s{x^y) around the two coordinate axes 
X, y with rotation angles a, (3 obtained from 

tana(x,2/) = -dys{x,y), tan/5(x,?/) = d:,s{x,y) (2) 

assuming that the rotation angles are small. Next, we translate the surface point along the 
z-axis from z = s{x^y) to z = 0. If the elevation is small, then we achieve this through a 
Taylor expansion of the stress tensor according to: 

Tij{x, y,z = s{x, y)) ^ Tij{x, y, 0) + s{x, y)d^Tij{x, y, z) (3) 

Whether the last equation is a good approximation also depends on the length of seismic 
waves (and therefore frequency) since Tij is the stress tensor associated with some displace- 
ment field. The greater the seismic wavelength in the vertical direction (along the z-axis), 
the smaller ^^r^j, resulting in a better approximation. Finally, analogous with the displace- 
ment field, the stress field will be divided into an unperturbed stress r^^ in case of a fiat 
surface and a stress tensor rf^ associated with the scattered field. Neglecting higher-order 
perturbations in s{x,y) and combining equations ([T]) and (|3]), we find 

^xz = -^(^. y)9zr^z + {d^s{x, y))T^^ + {dys{x, y))T^y (4) 
^yz = -^(^. y)dzr^z + {d^s{x, y))T^y + {dys{x, y))T^y (5) 
<z = -s{x^y)9zr^z (6) 

where all stress tensors being evaluated on the plane z = 0. We used r£(x, ^, 0) = because 
the unperturbed stress field produces a vanishing traction normal to a free, fiat surface. 
Linearizing the surface stress equations with respect to the surface profile s{x^y) is part of 
the Born approximation. 

We now have equations that allow us to calculate a stress field at z = in terms of a 
surface topography s{x^y) and unperturbed stress r^-. The two remaining problems are to 



calculate the unperturbed stress field r^j and its derivatives (at z = 0) given an unperturbed 
seismic field and then to calculate the scattered displacement field throughout the 
entire medium using boundary conditions that are derived from the effective surface loads 
rf^. For this purpose, it is convenient to express the unperturbed displacement field in terms 
of potentials 

^°(f , t) = V(/)(x, t) + V X if{x, t) (7) 

as we can now easily associate the scalar potential (j) with compressional modes and the 
vector potential with shear modes (V • ^0 = 0). For each of the potentials, the traction 
with respect to planes perpendicular to the z-axis reads 

CP =\^(i>8i, + 2iidA(t> (8) 

it" = 1^ (2a,(V x^)- A(e-; X ^))^ (9) 

where the subscript "i" in the last equation denotes the component i of the vector inside 
the brackets. A, ji are the Lame constants, and A is the Laplace operator. We will not give 
explicit expressions here for the potentials. They can be found in ^15^, J^. We just want to 
point out that in all cases studied in this paper, i. e. scattering of incident compressional, 
shear and fundamental Rayleigh waves, the scalar and vector potentials both need to be 
included due to partial mode conversion at refiection from the fiat surface (the Rayleigh 
field being a complex continuation of the bulk- wave refiection). 

The same equations ^ and ^ together with equation ^ and V - if = Q can be used 
to calculate the scattered waves. One starts with the surface load, equations ([i]) to ([g]), and 
solves for the amplitudes of the field potentials. In other words, the surface load forms a 
boundary condition for the scattered field. We can determine the potentials analytically only 
for very simple surface loads. In this paper, we calculate the amplitudes of plane waves con- 
tained in the scattered field as a function of the wave vector R that characterizes a sinusoidal 
surface profile and for different types of incident plane waves. Then the scattering problem 
is analogous to refiection of light from optical gratings, and results can be easily interpreted 
in wavenumber space. If the s{x^y) has a continuous spectrum, one needs to integrate over 
all possible horizontal wavenumbers of the scattered field, but because we consider load 
distributions rf^(x,^) with discrete spatial spectrum and since we understand the surface 
load as boundary condition of the scattered field (as opposed to a source distribution), no 
integration is required to calculate the amplitudes of the scattered waves. 
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As an example, let us consider an incident plane compressional wave. Its wave vector 
can be decomposed into a sum over the horizontal (along the plane z = 0) and vertical wave 
vectors: = fc^'^ + fc^'^. At reflection from a flat surface, the horizontal wave vector is 
preserved among all in and outgoing waves, but the vertical wave vector will have diflFerent 
values for the different wave modes. There is a simple sign flip in fc^'^ for the reflected 
compressional wave, but the shear wave generated at reflection has a vertical wavenumber 
^shr > k^'^ because the shear wave speed is smaller than the compressional wave speed. So 
in this case, the complete unperturbed displacement fleld including the incident and reflected 
waves is characterized by a single horizontal wave vector and three vertical wave vectors. 
According to equations ([s]) and ^ , the same wave vectors also characterize the unperturbed 
stress fleld. Substituting a sinusoidal function exp(±i^-x) for the surface topography ^(x, y) 
in equations Q to ([6]) , we flnd that the surface load is characterized by a discrete spectrum 
with horizontal wave vectors k^^ = fc°'^±^. It clearly follows that the scattered waves must 
have the same horizontal wave vectors. Therefore, in this simple example, the scattered fleld 
is composed of waves with four different wave vectors determined by the two k^^ and two 
different wave speeds a, /3 associated with compressional and shear waves. The scattered 
fleld can be represented in the form 

+S^(fc;'^)e^(^+(^)-^^-"^) + S^(fc!'^)e^(^^-(^>^^-"^) (10) 
where the 3D wave vectors k^ have vertical wavenumber 

kric) ^ ^^-{ktr (11) 

The scattering coefficients (7, S describe the scattering of the incident wave into compres- 
sional and shear waves. Small- wavelength sinusoidal surfaces, ^ scatter waves 
predominantly at large angles with respect to the specular reflection, whereas surfaces with 
^ ^ ^o,h gcg^i^i^gj. waves at small angles relative to the specular reflections. In general, the 
topographic spectrum will be continuous and extend to arbitrarily large wavenumbers 
Therefore, the wavenumber spectrum of the scattered fleld will extend to arbitrarily large 
horizontal wavenumbers fc^'^, and a bulk compressional or shear wave can scatter into a sur- 
face wave with imaginary vertical wavenumber and vice versa. Among the scattered surface 
waves, only the ones with a speciflc horizontal wavenumber correspond to the fundamental 
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Rayleigh mode, which means that scattering creates surface waves that are not Rayleigh 
waves. These waves cannot propagate freely because fundamental Rayleigh waves are the 
only surface waves that can propagate freely on homogenous media. In the context of New- 
tonian noise estimation, scattered surface waves that have wavelengths much shorter than 
Rayleigh waves can pose a problem. If these waves have significant amplitude, very dense 
seismic arrays would be necessary to monitor the surface displacement. 

The amplitudes of the scattered field at the origin x = are plotted in figures [l] and 
[2j Figure [T] also shows an example of a scattered displacement field with incident Rayleigh 
wave. In all cases, the amplitudes are calculated for incident waves with 1 m amplitudes 
(1 m vertical displacement for incident Rayleigh waves) and 1 m amplitude of the sinusoidal 
topography. The compressional-wave to shear-wave speed ratio is slightly smaller than 
2 for all plots in this section. Certain surface wavelengths Ag produce scattered waves 
with horizontal wavenumbers that correspond to the (real- valued) zeros of the Rayleigh 
function 

i?(fc^) = ((fc^)2 - {k^{f3)ff + 4(fc^)2fc^(a)fc^(/3) (12) 

The Rayleigh function appears in the denominator in coupling coefficients between seismic 
sources and fields near surfaces and so it also naturally emerges when solving the surface 
boundary equations ([S]) and ^ for the displacement fields. In general, the zeros of the 
Rayleigh function (usually called Rayleigh poles since R{k) is mostly found in denomina- 
tors) describe a resonance between compressional and shear waves near the surface. Near 
these points, the scattering coefficients go to infinity, and the Born approximation breaks 
down. This problem is partially a result of considering scattering from infinite sinusoidal 
topographies, which is certainly a rather artificial scenario. Often, the Rayleigh poles do not 
yield pathological behavior when integrating over a range of wavenumbers (see for example 
[19]). For the purpose of this paper, we will simply understand scattering into the Rayleigh 
poles as very efficient without further quantifying the scattering coefficients. Finally, we 
want to point out that the absolute value of the scattering coefficient is also a function of 
frequency. The plots in figures [T] and [2] were evaluated at / = 10 Hz. There is no simple fre- 
quency scaling valid for the entire wavenumber space, but scattering always becomes weaker 
at smaller frequencies if all other parameters are kept constant. 
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FIG. 1: The plot to the left shows an example of a scattered displacement field in the x,2:-plane. 
The incident wave is a Rayleigh wave with horizontal wavelength A-^^^^ ^ 500 m propagating along 
the x-direction. It is scattered from a sinusoidal surface that has a wavelength Ag = 27r//^ = 100 m. 
The incident wave has vertical amplitude 1 m at z = and the surface topography has amplitude 
1 m. In the Born approximation, the amplitude of the scattered wave scales linearly with the 
amplitude of the surface profile and linearly with the amplitude of the incident wave. One can see 
that the evanescent field is not of Rayleigh type because its phase is a function of z. The plot to 
the right shows the absolute value of the vertical displacement of the scattered field at the surface 
normalized by the amplitudes of the topography and incident wave. The topographic wave vector 
is expressed in units of the horizontal wavenumber of the incident Rayleigh wave. A resonance 
effect between scattered shear and compressional waves causes the scattering amplitude to go to 
infinity for certain which lie on the dark colored circle. The Born approximation breaks down 
at these wave vectors, which are known as the Rayleigh poles. For this reason, it is not possible 
to calculate a total integrated scatter for harmonic topographies in the Born approximation. The 
problem seems to be analogous to the emergence of secular terms in perturbation theory that is 
often solved by including dispersion terms in the solution. 

III. TOPOGRAPHIC DATA AND SITE SELECTION 



As noted, we use information about surface topography to investigate seismic scattering at 
two representative sites. As the arms of current GW detectors are on the order of kilometers, 
we will analyze locations of 10km x 10km area. We have acquired surface topography data 
covering the United States from the USGS National Elevation Dataset. The elevation data 
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FIG. 2: The plots show absolute values of the scattering coefficients into vertical surface displace- 
ment for incident compressional (left) and shear (right) waves propagating along the x-direction 
as a function of wave vector /? of a sinusoidal surface. The topographic wave vector is expressed 
in units of the horizontal wavenumber of the incident compressional or shear wave. The incident 
shear wave is in so-called SV polarization, which has no displacement component along the hori- 
zontal direction. In contrast to Rayleigh-wave scattering as shown in figure [TJ bulk- wave scattering 
coefficients also exhibit forward scattering into the Rayleigh poles that lie on the dark circle. 

is provided in blocks of 1° latitude by 1° longitude with 1 second resolution. This high 
resolution data is necessary to understand how rapidly the elevation changes over the course 
of land data on the scales of which we are concerned in the context of topographic scattering. 

As blocks of latitude and longitude are neither regularly spaced nor equal area, it is 
necessary to perform a coordinate transition to analyze blocks of the same area. Therefore, to 
map the elevation data, we use Lambert's equal-area cylindrical projection. After projection, 
the data is interpolated to be on a regular 10km x 10km grid covering the country. The 
benefits of this grid are that the arms of current GW detectors are on the order of kilometers, 
and so it is fair to think of each block in the grid as a location for a future GW detector. 
We provide two plots in figure |3) The top one is the mean elevation at each grid point. The 
bottom plot shows the root-mean-square (rms) of the elevation at each grid point. This map 
will be used to identify locations of low rms. This is important as locations with rapidly 
changing elevation will show stronger scattering. Therefore, we desire to pick locations low 
in elevation rms. 

However, topography is not the only determinant of the complexity of ambient seismic 



9 



^0 1000 2000 3000 4000 5000 

X coordinate [km] 




1000 2000 3000 4000 5000 

X coordinate [km] 

FIG. 3: Tlie plot on tlie top sliows tlie elevation map of the US. Although not further studied 
in this paper, one important observation is that seismic noise especially above 1 Hz is weaker in 
regions with high elevation. In terms of site selection, the challenge is to find a high-elevation site 
(i. e. low seismic noise) with small changes in elevation over an area of about lOkmxlOkm to 
minimize scattering. For this purpose, the bottom map shows the rms of the elevation calculated 
for areas of this size. Every point of the plot represents one of these squares so any light colors 
in the rms plot within high-elevation areas is a possible candidate for a low-seismic noise, low 
scattering site. 

fields, and this must also be factored when choosing sites. In general, flat land is better 
populated than land that has rapidly changing elevation. Therefore, land outside of the 
mountains has the shortcoming of increased anthropogenic seismic noise. For this reason, 
we desire a location that has low-rms and is therefore flat, while also seismically quiet. We 
use the rms map as a flrst flgure-of-merit to identify promising sites, and these locations 
will be examined more carefully with detailed scattering calculations in the next section. 
Their elevation maps are shown in flgure [4} One site represents a high-rms area, the other 
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FIG. 4: The two maps display elevation data of low-rms and high-rms locations. A low-rms 
location with weak seismic noise would be an ideal site candidate for a future GW detector. Here, 
we only study the scattering problem, but both sites also have known seismic noise measured by 
the USArray stations F13A for the map to the left, and K07A for the map to the right. Both sites 
have weak seismic noise at frequencies between lOmHz and 10 Hz. 

one a low-rms area. The maps are centered within a few km to seismic stations of the 
USArray to ensure that we have local seismic data. For example, the low-rms site is near 
the station K07A of the transferrable array (TA), where noise spectra close to the low-noise 
model are observed (see [20j for the definition of the seismic-noise models). The center of 
its map shown in figure [4] is located at N42.64, W119.13 in Oregon close to the border with 
California. The location of the station F13A coincides with the center of the high-rms site 
at N45.79, W114.33 in Montana. It also represents a low-noise site. As can be seen from 
the map, K07A has a peak-to-peak elevation change of about 14 m over 5 km distance to the 
center in any direction, which is much less than we had initially anticipated for any high- 
elevation site. Together with the seismic spectrum, as shown in figure [5| it is demonstrated 
that low-noise, fiat surface sites exist in the US that could be considered as locations for 
a future GW detector. We want to emphasize though that there are many more sites that 
can be found with similar characteristics. The Oregon site was just the best one according 
to our simplified rms criterion. In addition, it is certainly true that topographic scattering 
and seismic noise are only two of a large number of site-selection criteria. 
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FIG. 5: The plot shows the seismic spectral histogram measured at station K07A of the trans- 
ferrable array (part of the USArray) based on two years of 3-hour averaged spectra. The two 
dashed lines are the high-noise and low- noise global seismic models. A more careful seismic anal- 
ysis needs to be carried out to fully characterize a site, but the spectrum proves that flat surface, 
low- noise sites exist in the US. 

IV. PHENOMENOLOGICAL SCATTERING ANALYSIS FOR F13A AND K07A 

The last step is to expand the topographic maps of the two selected sites into sums 
over sinusoidal surface perturbations and to calculate the amplitudes of scattered waves. 
Because the mode composition of the seismic fields at the sites is unknown, we assume that 
the unperturbed seismic field is composed of fundamental Rayleigh waves. This simplified 
model is sufficient to discuss how topographic scattering affects Newtonian-noise estimation 
and subtraction in future GW detectors. 

Figure |6] shows the histograms of ID spectra of the topographic maps along all grid lines in 
two orthogonal directions. As scattering coefficients in wavenumber space are proportional 
to the topographic spectrum, the widths of the spectral histograms give an idea of how much 
variation in scattering is to be expected as a function of propagation direction of the incident 
seismic wave. The units of the spectra are chosen such that a simple sum over the harmonic 
components yield the ID topography along grid lines. In this way, we can test easily if the 
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FIG. 6: The spectral histograms are derived from ID Fourier transforms of the topographic maps 
(figure [4]) along all grid lines in two orthogonal horizontal directions. The white dashed curves are 
the 90th percentiles of the histogram. 

conditions for the Born approximation still hold. For example, scattering of a 10 m long 
Rayleigh wave from 1 km long topographic features at site F13A cannot be described in the 
Born approximation since the corresponding topographic amplitude is larger than the length 
of the incident wave. Comparing the two histograms we find that the amplitudes at larger 
wavenumbers have similar values whereas the amplitudes at small wavenumbers diflFer by 
almost two orders of magnitude. In general, we should always expect topographic spectra to 
decrease with increasing n (the higher the mountain the wider its base). Therefore, small- 
scattering should be more pronounced. 

Next, we multiply the scattering coefficient shown in figure [T] with the 2D topographic 
spectra of each site. The result is the amplitude of scattered waves as a function of topo- 
graphic wavenumber normalized by the amplitude of the incident Rayleigh wave. For this 
analysis, we chose to calculate the scattering at lower frequency, i. e. at 3 Hz, and to be 
more careful about the estimate of compressional and shear- wave speeds a, /3 at the two 
sites (which also determine the Rayleigh- wave speed). According to USGS geologic maps 
( http: / / tin.er.usgs.gov / geology / state /[ ) , the geologic unit at the low-rms Oregon site K07A 
is "unconsolidated to semi-consolidated clay, silt, sand or gravel". Speed of compressional 
waves in these materials is typically a = 1800 m/s with comparatively low shear- wave speeds. 
We assume /3 = 900 m/s, which falls into the shear- wave range specified in [21j for stiff soil. 
F13A lies in the Idaho Batholith, which contains "faintly gneissic quartz monzonite, gran- 
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odiorite, and similar rocks" . Depending on the state of weathering of the rock, compressional 
wave speed can range from 2000 m/s to 5000 m/s. Because the area also contains alluvial 
deposits, we chose a = 2500 m/s near the lower end of the range and shear- wave speed 
/5 = 1800 m/s. The scattering coefficients are shown in figure [?[ As expected, scattering is 
significantly stronger at F13A. As the (horizontal) wave vector of a scattered wave is simply 
the sum of the wave vector of the incident Rayleigh wave propagating along the x-direction 
and the topographic wave vector ^, the plots can be understood as the spatial spectra of 
a seismic wave that entered the site from the west as a plane wave and has just left it in 
the east. As we mentioned before, scattering of Rayleigh waves into Rayleigh waves has to 
be much weaker than suggested by our results since the Born approximation breaks down 
near the Rayleigh poles giving rise to infinite scattering coefficients. Neglecting scattering 
at these wave vectors, the main effect of irregular surface topography in the Born approxi- 
mation is to broaden narrow features in the spatial spectrum of the incident wave. Due to 
the Rayleigh poles, this spreading should occur into arc like shapes. 

We want to point out that scattering into Rayleigh waves does not inffuence the perfor- 
mance of a NN subtraction filter. The seismic array is already designed to allow subtraction 
of NN from Rayleigh waves. It does not make any difference whether these waves are in- 
cident from outside or produced by scattering. However, broadening of narrow features in 
wavenumber spectra can potentially limit NN subtraction because it would be challenging 
to design a seismic array that provides data to correctly estimate seismic amplitudes over 
continuous intervals of wavenumbers (at each frequency). Even though a total integrated 
scatter cannot be obtained due to the Rayleigh poles, scatter coefficients can be summed if 
one excludes a ring that lies on the poles. Excluding a ring of width dk = 0.01 /cr, we find 
a total integrated scatter of 0.04 for the F13A site. Similarly, the total integrated scatter at 
the low-rms site K07A is 0.002. If the goal is to reduce NN by a factor 100, then scattering 
can not be neglected at the F13A site. 

V. CONCLUSION 

We discussed topographic scattering in the context of GW-detector site selection and 
Newtonian-noise subtraction. We found that the total contribution of waves scattered from 
topography can exceed values of 0.01, which makes topographic scattering relevant to NN 
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FIG. 7: The contour plots show the scattering coefficients for incident Rayleigh waves. Very small 
surface wavenumbers cause a weak broadening of seismic spatial spectra. Scattering from very 
high surface wavenumbers can potentially have large effects on spatial spectra of seismic fields, 
but the corresponding scattered waves tend to have negligible amplitudes as can be seen in the 
two plots (even at the high-rms site). The wavenumber range between these two extremes is 
most interesting, and it also contains the poles of the scattering coefficients in the case of incident 
Rayleigh waves. Even though a total integrated scatter cannot be obtained due to the Rayleigh 
poles, scatter coefficients can be summed if one excludes a ring that lies on the poles. The results 
obtained for a ring of width d/c = 0.01 /cr are 0.002 for the low-rms site and 0.04 for the high-rms 
site. 

subtraction in future low-frequency GW detectors assuming NN reduction by a factor 100. 
Although scattering into Rayleigh waves cannot be described in the Born approximation, 
we explained that this relatively efficient scattering channel would not pose a problem for 
NN subtraction independent of the value of the scattering coefficient. Since it is possible 
that scattering into Rayleigh waves significantly dissipates energy from the incident wave at 
strong scattering sites, the absolute values of scattering coefficients may be overestimated 
though. This issue can for example be resolved with numerical simulations of scattering into 
Rayleigh waves. For scattering into wavenumbers that do not correspond to Rayleigh poles, 
the outlined theory provides robust predictions of scattering coefficients. Since scattering 
into body fields is negligible for all types of incident waves, the problem of topographic 
scattering can be avoided by constructing the detector sufficiently deep underground. This 
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scenario would have to be studied carefully though since scattered surface waves are not 
generally of Rayleigh type and can therefore generate displacement deeper than known for 
Rayleigh fields. 

The main challenge with subtraction of NN from surface waves is the design of the 
seismic array. The requirement for the array diameter is derived from the longest seismic 
wavelengths in the field at lowest frequency, whereas the array density, which will be variable 
throughout the array, depends on the shortest wavelengths at highest frequency. In fields 
where there is no unique correspondence between frequency and wavelength, it will be very 
difficult to design a seismic array for NN subtraction. As was shown, this situation occurs 
when a site exhibits strong topographic scattering that causes amplitudes in wavenumber 
space to spread over adjacent wave vectors. This spreading happens with relatively high 
scattering amplitudes because it is associated with the lowest wavenumber components of 
the topography that typically have the largest amplitudes. We have not repeated the same 
analysis for scattering of bulk waves since the effects are very similar. In addition, scattering 
of bulk waves into Rayleigh poles is generally driven by higher- wavenumber components of 
the surface topography that have weaker amplitudes. 

In the future, a simulation of topographic scattering should be carried out that includes 
the seismometer array and the noise subtraction. Then one can address questions such as 
how the array needs to be adjusted to minimize the effect of small- wavelength scattering on 
subtraction residuals. Although we found that topographic scattering can be negligible at 
fiat sites or that it does not lead to adverse effects, it is certainly desirable to understand 
how to design arrays to achieve subtraction goals at arbitrary sites. 
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